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O- , Abstract 

We describe the program HFBTHO for axially deformed configurational Hartree-Fock- 
Q \ Bogoliubov calculations with Skyrme-forces and zero-range pairing interaction using Harmonic- 
Oscillator and/or Transformed Harmonic-Oscillator states. The particle-number symmetry is 
approximately restored using the Lipkin-Nogami prescription, followed by an exact particle 
number projection after the variation. The program can be used in a variety of applications, 
including systematic studies of wide ranges of nuclei, both spherical and axially deformed, 
^ ' extending all the way out to nucleon drip lines. 
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>• . Program Summary 

5h ' Title of the program: HFBTHO (vl.66p) 
d | 

Catalogue number: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Program summary URL: 

Licensing provisions: none 

Computers on which the program has been tested: Pentium-Ill, Pentium-IV, AMD-Athlon, IBM 
Power 3, IBM Power 4, Intel Xeon 

Operating systems: LINUX, Windows 

Programming language used: FORTRAN-95 
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Memory required to execute with typical data: 59 MB when using N s h = 20 

No. of bits in a word: 64 

No. of processors used: 1 

Has the code been vectorized?: No 

No. of bytes in distributed program, including test data, etc.: 

No. of lines in distributed program: 7876 lines 

Nature of physical problem: Solution of self-consistent mean-field equations for weakly bound 
paired nuclei requires correct description of asymptotic properties of nuclear quasiparticle wave 
functions. In the present implementation, this is achieved by using the single-particle wave 
functions of the Transformed Harmonic Oscillator, which allows for an accurate description of 
deformation effects and pairing correlations in nuclei arbitrarily close to the particle drip lines. 

Method of solution: The program uses axial Transformed Harmonic Oscillator single-particle 
basis to expand quasiparticle wave functions. It iteratively diagonalizes the Hartree-Fock- 
Bogolyubov Hamiltonian based on the Skyrme forces and zero-range pairing interaction until 
the self-consistent solution is achieved. 

Restrictions on the complexity of the problem: Axial-, time-reversal-, and space-inversion sym- 
metries are assumed. Only quasiparticle vacua of even-even nuclei can be calculated. 

Typical running time: 4 seconds per iteration on an Intel Xeon 2.8 GHz processor when using 
N sh = 20 

Unusual features of the program: none 
PACS: 07.05.T, 21.60.-n, 21.60.Jz 
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1 Introduction 

Nuclear structure theory strives to build a comprehensive microscopic framework in which 
bulk nuclear properties, nuclear excitations, and nuclear reactions can all be described. Exotic 
radioactive nuclei are the critical new focus in this quest. The extreme isospin of these nuclei 
and their weak binding bring new phenomena that amplify important features of the nuclear 
many-body problem. 

A proper theoretical description of such weakly bound systems requires a careful treatment 
of the asymptotic part of the nucleonic density. An appropriate framework for these calculations 



is Hartree-Fock-Bogoliubov (HFB) theory, solved in coordinate representation PQ |2]. This 
method has been used extensively in the treatment of spherical nuclei [3] but is much more 
difficult to implement for systems with deformed equilibrium shapes. There have been three 
ways of implementing deformation effects into the coordinate-space HFB. The oldest method, 
the so-called two-basis method jUEHH], is based on the diagonalization of the particle-particle 
part of the HFB Hamiltonian in the self-consistent basis, obtained by solving the HF problem 
with box boundary conditions. The disadvantage of this method is the appearance of a large 
number of positive-energy free-particle (box) states, which limits the number of discretized 
continuum states (the maximum single-particle energy taken in this method is usually less 
than lOMeV). 

The second, very promising strategy, the so-called canonical-basis HFB method, utilizes the 
spatially localized eigenstates of the one-body density matrix without explicitly going to the 
quasiparticle representation 0|H1IH|. Finally, an approach to axial coordinate-space HFB has 
recently been developed that uses a basis-spline method [TUHH]. While precise, these two latter 
methods are not easy to implement and, because they are time-consuming, cannot be used in 
large-scale calculations in which a crucial factor is the ability to perform quick calculations for 
many nuclei. 

In the absence of fast coordinate-space solutions to the deformed HFB equations, it is useful 
to consider instead the configuration-space approach, whereby the HFB solution is expanded in 
some single-particle basis. In this context, the basis of a harmonic oscillator (HO) turned out 
to be particularly useful. Over the years, many configuration-space HFB+HO codes have been 
developed, either employing Skyrme forces or the Gogny effective interaction [T2 l IT8 | IT U 1X5 } IT6] . 
or using a relativistic Lagrangian ^2j in the context of the relativistic Hartree-Bogoliubov 
theory. For nuclei at the drip lines, however, the HFB+HO expansion converges slowly as a 
function of the number of oscillator shells [3], producing wave functions that decay too rapidly 
at large distances. 

A related alternative approach that has recently been proposed is to expand the quasipar- 
ticle HFB wave functions in a complete set of transformed harmonic oscillator (THO) basis 
states [IB], obtained by applying a local-scaling coordinate transformation (LST) [THJ |20] to 
the standard HO basis. Applications of this HFB+THO methodology have been reported both 
in the non-relativistic [21J and relativistic domains [22J. In all of these calculations, specific 
global parameterizations were employed for the scalar LST function that defines the THO ba- 
sis. There are several limitations in such an approach, however. For example, the minimization 
procedure that is needed in such an approach to optimally define the basis parameters is com- 
putationally very time-consuming, making it very difficult to apply the method systematically 
to nuclei across the periodic table. 

Recently, a new prescription for choosing the THO basis has been proposed and employed in 
self-consistent large-scale calculations [231- For a given nucleus, the new prescription requires as 
input the results from a relatively simple HFB+HO calculation, with no variational optimiza- 
tion. The resulting THO basis leads to HFB+THO results that almost exactly reproduce the 
coordinate-space HFB results for spherical nuclei [21]. Because the new prescription requires 
no variational optimization of the LST function, it can be applied in systematic studies of 
nuclear properties. In order to correct for the particle number nonconservation inherent to the 
HFB approach, the Lipkin-Nogami prescription for an approximate particle number projection, 
followed by an exact particle number projection after the variation has been implemented the 



code HFBTHO (vl.66p) 

The paper is organized as follows. Section El gives a brief summary of the HFB formalism. 
The implementation of the method to the case of the Skyrme energy density functional is 
discussed in Sec. H3 together with the overview of the THO method and the treatment of 
pairing. Section 0] describes the code HFBTHO (vl.66p). Finally, conclusions are given in 
Sec. El 

2 Hartree-Fock-Bogoliubov Method 

A two-body Hamiltonian of a system of fermions can be expressed in terms of a set of annihi- 
lation and creation operators (c,^): 

ti = / J e n 1 n 2 C m C n2 ' 4 / , Vnin 2 n 3 n4, C ni C n 2 C «4 C n.3' {*■) 

niri2 nin2n3n4 

where v nin2ri3n4 = (n 1 n 2 |V A |n 3 n4 — n^n 3 ) are anti-symmetrized two-body interaction matrix- 
elements. In the HFB method, the ground-state wave function |$) is defined as the quasiparticle 
vacuum «fe|$) = 0, where the quasiparticle operators (a,a^) are connected to the original 
particle operators via the linear Bogoliubov transformation 
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which can be rewritten in the matrix form as 

5)-(££)U)- m 

Matrices U and V satisfy the relations: 

rfu + v*v = i, UU ] + V*V T = 1, U T V + V T U = 0, UV* + V*U T = 0. (4) 
In terms of the normal p and pairing k one-body density matrices, defined as 

Pnn , = ($|ci,c n |$) = (V*V T ) nn ,, n nn , = ($|c n /c n |$) = (V*U T ) nn ,, (5) 

the expectation value of the Hamiltonian (JTJ) is expressed as an energy functional 

E[p > K] = T^T = Tr [(e + |r)p] " |Tr [Ak * ] ' (6) 

where 

-*- nin3 / ^niri2n3n4Pn4ri2 5 ^niri2 2 / ., ^'nin2n3n4"'n3n4 - V / 

ri2n4 «3"4 

Variation of energy (JHJ) with respect to p and k results in the HFB equations: 

e+T-A A \fU\ fU 

-A* -(e + r)* + A J \ V J \V 

where the Lagrange multiplier A has been introduced to fix the correct average particle number. 
It should be stressed that the modern energy functionals (JHJ) contain terms that cannot be 
simply related to some prescribed effective interaction, see e.g., Ref. J2Z1I2E! for details. In this 
respect the functional (jHJ) should be considered in the broader context of the energy density 
functional theory. 



3 Skyrme Hartree-Fock-Bogoliubov Method 

3.1 Skyrme Energy Density Functional 

For Skyrme forces, the HFB energy © has the form of a local energy density functional, 

E[p,p] = JdhH(r), (9) 

where 

H{r) = H(r) + H(r) (10) 

is the sum of the mean-field and pairing energy densities. In the present implementation, we 
use the following explicit forms: 

H(r) = fer + i* [ {l + \x Q )p*-{\ + x Q )Y. P 2 q ] 

+ ¥i [ (1 + §zi) p (r - | Ap)] - (| + xi) E P, (r« - f Ap,)] 
+ |*2 [ (1 + |x 2 ) p (r + ±Ap) - (| + x 2 ) E P 9 {r q + \A Pq )} 

+ ^3P^(l + ix 3 )p 2 -(x3 + |)E^ " (U) 



I (tixi + t 2 x 2 ) e j?,- + 1 (*i - * 2 ) e j;,« 

|Wb E e iiA [pVfeJij + E PgVfcJg,jj] 



and 

" 7 



ff(r) = 1V„ 



Po 



£# ( 12 ) 
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Index q labels the neutron (q = n) or proton (q = p) densities, while densities without index q 
denote the sums of proton and neutron densities. H(r) and H(r) depend on the particle local 
density p(r), pairing local density p(r), kinetic energy density r(r), and spin-current density 
Jij(r): 

p(r) = p(r,r), p(r) = p(r,r), 

(13) 
r(r) = V r V r ,p(r,r')U r , J y (r) = £ (V, - V$) /^(r, rQ Ur , 

where p(r, r'), pi(r, r'), p(r, r'), pj(r,r') are defined by the spin-dependent one-body density 
matrices in the standard way: 

p(rcr, rV) = ~p(r, r')<W + § E(o-|o-i|o-')Pi( r > r ')> 

(14) 
p(ra, rV) = |p(r, r')S aa > + \ '^(a\ai\a')p i (r, r'). 

i 

We use the pairing density matrix p, 

p(r<7, rV) = — 2cr'«;(r,cr, r', — a'), (15) 



instead of the pairing tensor k. This is convenient when describing time-even quasiparticle 
states when both p and p are hermitian and time-even |2]. In the pairing energy density (|12|). 
we have restricted our consideration to contact delta pairing forces in order to reduce the 
complexity of the general expressions 0I2EJ- 

3.2 Skyrme Hartree-Fock-Bogoliubov Equations 

Variation of the energy Q with respect to p and p results in the Skyrme HFB equations: 
h(r,a,a') h{r,cr,cr') \ ( U{E, rcr') \ ( E + X \ { U(E,ra 



a , \h(r,a,a') -h(r,a,a>) J\ V(E,ra') J \ E - X J \ V(E, ra) J ' l ° J 

where local fields h(r, a, a') and h(r, a, a') can be easily calculated in the coordinate space by 
using the following explicit expressions: 

h q {r, a, a') = -VM q V + U q + ^ (V^^- + B^-V^) , 
h q (r,a,a') = V (l-vJ£\ \ p q , 
where 

M 1 = <£ + 1*1 [(1 + 5*0 P ~ (*1 + I) P?] + 1*2 [(1 + |^2) P + (X 2 + |) pj] , 

£g,ij = ~\ (tlXi +t 2 ^2) Jy + \{t\ ~h) Jq,ij + ^W J2 £ ijkV k (p + Pq,) , 

ijk 

U q = t [(l + 2^0) P ~ {X + 5) Pq] 

+ \t x [ (1 + \x x ) (t - |Ap) - (an + §) (r, - §ApJ] 

(18) 

+ \h [ (1 + \x 2 ) (r + §Ap) + (x 2 + \) (r, + ±Ap g )] 

+ ^3P a [(l + ^ 3 )(2 + «)p-2(x 3 + i)p (? +(l-X3)f]-^(^) 7 EP 



~,2 

g 
q 



± (tlXl + t 2 X 2 ) J2 1% + \{tl~ t 2 ) £ Jg,y - I Wo £ %fcV fc [Jy + J g ,y] . 
ij q,ij ijk 

Properties of the HFB equation in the spatial coordinates, Eq. (jTUjl . have been discussed in 
Ref . [2] . In particular, it has been shown that the spectrum of eigenenergies E is continuous for 
|J5|>— A and discrete for \E\<—X. In the present implementation, we solve the HFB equations 
by expanding quasiparticle wave functions on a finite basis; therefore, the quasiparticle spectrum 
Ek becomes discretized. Hence in the following we use the notation V^(ra) = V(E k ,ra) and 
CTjfc(rcr) = U(E k ,ra). Since for E k >0 and A<0 the lower components Vfc(rcr) are localized 
functions of r, the density matrices, 

p(ra,rV) = £V*(r«7)Vjr(rV), (19) 

k 

p(ra,rV) = - J>,(rcr)?7 fe *(rV'), (20) 

k 
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are always localized. The orthogonality relation for the single-quasiparticle HFB wave functions 
reads 

[d 3 rY^[Uk(™)U k ,(ra)+V£(ra)V k ,(ra)]=5 k>k/} (21) 

and the norms of lower components N k , 

N k = /"d 3 r^|V fe (ra)| 2 , (22) 

define the total number of particles 

N= fd 3 r p(r) = Y,N k . (23) 

3.3 Axially Deformed Nuclei 

For spherical nuclei, Skyrme HFB equations are best solved in the coordinate space, because 
Eq. ([TE]) reduces in this case to a set of radial differential equations [22! • I n the case of deformed 
nuclei, however, the solution of a deformed HFB equation in coordinate space is a difficult and 
time-consuming task. For this reason, here we use the method proposed by Vautherin [30 , 
which combines two different representations. The solution of the deformed HFB equation is 
carried out by diagonalizing the HFB hamiltonian in the configurational space of wave- functions 
with appropriate symmetry, while evaluation of the potentials and densities is performed in 
the coordinate space. Such a method is applicable to nonaxial deformations [IE] , but typical 
computation time for large-scale mass-table calculations is prohibitively large. In the present 
implementation, we make the restriction to axially-symmetric and reflection-symmetric shapes 
in order to obtain HFB solutions within a much shorter CPU time. 

In the case of axial symmetry, the third component J z of the total angular momentum is 
conserved and provides a good quantum number Q k . Therefore, quasiparticle HFB states can 
be written in the following form: 



U k (r,a,r) . 
V k (r,<r,T) I wl,) 



U k (r,z) \ A - l(n \i( U k (r,z) \ _*a+ Vv , - 

V k + (r, z) ) 6 X+1 ' 2[a) + I V k ~(r, z^ J e x - 1 / 2 ^ 



, (24) 



where A 1 * 1 = Q k ± 1/2 and r, z, and if are the standard cylindrical coordinates defining the 
three-dimensional position vector as r = (r cos ip,r sin ip,z), while z is the chosen symmetry 
axis. The quasiparticle states ()24|) are also assumed to be eigenstates of the third component 
of the isospin operator with eigenvalues q k = +1/2 for protons and q k = —1/2 for neutrons. 

By substituting ansatz ([23]) m t° Eq. ([TE]) . the HFB equation reduces to a system of equations 
involving the cylindrical variables r and z only. The same is also true for the local densities, 



i.e., 



p(r,z) -- 


= E 

k 


r(r,z) -- 


- E 



(l^ + M)l 2 + l^-M)| 2 ), 

+ |V 2 V; + (r,z)| 2 + |V 2 \4-(r^)| 2 + l|A-\4 + (r,,)| 2 ^ 
(V • J)(r, z) = J] (v r V^(r, z)V z V k {r, z) + ^V k + {r, z) [V r V+(r, z) - V z V k {r, z)] 

-V r V k -(r,z)V z V k + (r,z) - ^V k {r )Z ) [V r V k {r ) z) + V z V k + {r ) z)\^ 
~p(r,z) = -J2 {V k + (r,z)U+(r,z) + V k -(r,z)U k (r,z)), 

k 

(25) 
where V r = d/dr and V 2 = d/dz. In addition, when tensor forces are considered, the following 
additional densities have to be calculated: 

J r ^z) = J2 {V r V k + (r,z)V k -(r,z)-V r V k -(r,z)V k + (r,z)), 

k 

J<pr(r,z) = J2 (—V k + (r,z)V k -(r,z) + —V k -(r,z)V k + ( ri z) 
k \ r r 

J zv (r,z) = Y, (VzV k + (r,z)V k -(r,z)-V z V k -(r,z)V k + (r,z)). 

k 

J^(r,z) = J2 (—V k + (r,z)V k + (r,z)-—V k -(r,z)V k -(r,z) 



(26) 



k 

where indices denote the cylindrical components of the tensor J^, while all remaining compo- 
nents vanish due to the cylindrical symmetry, i.e., J rr (r, z) = J zz {r, z) = J w (r, z) = J rz (r, z) = 
J zr (r,z)=Q 

Due to the time- reversal symmetry, if the fcth state, defined by the set {U k , Uj~, V k , V k ~, fik}, 
satisfies the HFB equation (|TE|). then the kth state, corresponding to the set defined by 
{U k , —U k , V k + , —V k ~, — fik}, also satisfies the HFB equation for the same quasiparticle energy 
E k . Moreover, all wave functions in cylindrical coordinates are real. Contributions of time- 
reversal states k and k are identical (we assume that the set of occupied states is invariant with 
respect to the time-reversal), and we can restrict all summations to positive values of fik while 
multiplying total results by a factor two. In a similar way, one can see that due to the assumed 
reflection symmetry, only positive values of z need to be considered. 

3.4 HO and THO Wave Functions 

The solution of the HFB equation (J16J1 is obtained by expanding the quasiparticle function (|24j) 
in a given complete set of basis wave functions that conserve axial symmetry and parity. The 
program HFBTHO (vl.66p) is able to do so for the two basis sets of wave functions: HO and 
THO. 



The HO set consists of eigenf unctions of a single-particle Hamiltonian for an axially de- 
formed harmonic oscillator potential. By using the standard oscillator constants: 

and auxiliary variables 

£ = zp z , V = r 2 /3^ (28) 

the HO eigenfunctions are written explicitly as 

$ Q (r,a) = ^ r ( r )^ nz ( z ) X x(a), (29) 

V Z7T 



(30) 



where 



^n.{z) = Pl'^nM = N n jl /2 e-?/ 2 H nM (Z). 

Hn z {£.) and L^ r {rj) denote the Hermite and associated Laguerre polynomials [SI], respectively, 
and the normalization factors read 

( 1 \ 1/2 » / n l \ 1/2 
JV„ Z = and A^ A = ^^ . (31) 

™ 2 Vv^2^nj; n " VK+|A|)7 

The set of quantum numbers a = {n r ,n z ,A, £} includes the numbers of nodes, n r and n z , in 
the r and 2 directions, respectively, and the projections on the z axis, A and E, of the angular 
momentum operator and the spin. 

The HO energy associated with the HO state (|25|l reads 



e a = (2n r + |A| + l)tko± + (n z + \)fvuj z , (32) 

and the basis used by the code consists of M =(N sh +l)(N sh +2)(N sh +3)/6 states having the 
lowest energies e a for the given frequencies fkv± and ftw z . In this way, for the spherical basis, 
i.e., for ftw±=huj z , all HO shells with the numbers of quanta A^=0. . . N s h are included in the 
basis. When the basis becomes deformed, fkj±^hw z , the code selects the lowest-HO-energy 
basis states by checking the HO energies of all states up to 50 HO quanta. Note that in this 
case the maximum value of the quantum number Q^, and the number of blocks in which the 
HFB equation is diagonalized, see Sec. 13.51 depend on the deformation of the basis. 

The THO set of basis wave functions consists of transformed harmonic oscillator functions, 
which are generated by applying the local scale transformation (LST) [P9*l E21 120] to the HO 
single-particle wave functions ()29|). In the axially deformed case, the LST acts only on the 
cylindrical coordinates r and z, i.e., 



r'(r, z) = r ^— , 
z'(r,z) =z*yp-, 



(33) 



and the resulting THO wave functions read 



*«<-• *> - tWW< (^) *- &w) £**<*>• < 34 > 




where 

*=f M + W (35) 

and f(7Z) is a scalar LST function. In the code HFBTHO (vl.66p), function f(TZ) is chosen as 
in Ref. [23J. It transforms the incorrect Gaussian asymptotic behavior of deformed HO wave 
functions into the correct exponential form. Below, we keep the same notation $o,(r, a) for both 
HO and THO wave functions, because expressions in which they enter are almost identical in 
both cases and are valid for both HO and THO variants. 

3.5 HFB Diagonalization in Configurational Space 



We use the same basis wave functions to expand upper and lower components of the quasipar- 

U k {r,a,r) = Xq k {^)^2U ka ^> a {r,a), 



tide states, i.e., 



Vfc(r,a,r) = Xq k ( T )^2 v ka^a(r,a) } 



(36) 



where $ a (r, a) are the HO or THO basis states. Note that the same basis $ a (r, a) is used for 

protons and neutrons. 

Inserting expression (jHH|l into the HFB equation (fTT)|) and using the orthogonality of the basis 

states, we find that the expansion coefficients have to be eigenvectors of the HFB Hamiltonian 

matrix 

fcta*) _ \(ik) hiiu) \ ( u k \ _ / U k 



~ h M _/,fe) + A (») )\V k )- Ek \V k ) ' (37) 

where the quasiparticle energies E k , the chemical potential \( qk >, and the matrices 

h% = ($aN<^> and h® = ($a\h q \$f3) (38) 

are defined for a given proton (q k =+l/2) or neutron (q k —— 1/2) block. 

Proton and neutron blocks are decoupled and can be diagonalized separately. Furthermore, 
in the case of axially deformed nuclei considered here, fl k =A k +T lk is a good quantum number 
and, therefore, matrices h£p and h£l are block diagonal, each block being characterized by a 
given value of Q. Moreover, for the case of conserved parity considered here, 7r=(— 1)™ Z+A is 
also a good quantum number, and each of the Q k blocks falls into two sub-blocks characterized 
by the values of 7r=±l. Finally, due to the time-reversal symmetry, the Hamiltonian matrices 
need to be constructed for positive values of Q k only. 

3.6 Calculations of Matrix Elements 

As discussed in Sect. 13.21 local densities (j25|l and average fields, (|T7j) and (|TH)l . are calculated 
in the coordinate space. Therefore, calculation of matrix elements (138)) amounts to calculating 
appropriate spatial integrals in the cylindrical coordinates r and z. In practice, the integration 
is carried out by using the Gauss quadratures J2I] for 22 Gauss-Hermite points in the z > 
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direction and 22 Gauss-Laguerre points in the r direction. This gives a sufficient accuracy for 
calculations up to N sh = 40. 

In the case of the HO basis functions, the integration is performed by using the Gauss 
integration points, £ n and ri m , for which the local densities and fields have to be calculated at 
the mesh points of z n =b z ^ n and r m =b±r]rn ■ As an example, consider the following diagonal 
matrix element of the potential U q (r, z) (fTSjl. 



U q 

w act 



A 2, 



dz / rdr U q (z,r) ^J nz (z)^J nr (r 



(39) 



Inserting here the HO functions if>n x ( z ) an< ^ ^it ( r ) (EHJ), and changing the integration variables 
to dimensionless variables £ and rj, the above matrix element reads 



U q 







(40) 



where 

u q fov) = *v q (Zb x ,Snb ± ). (4i) 

Here, the Gauss integration quadratures can be directly applied, because the HO wave functions 
contain appropriate exponential profile functions. 

The situation is a little bit more complicated in the case of the THO basis states where, 
before calculating, one has to change variables with respect to the LST functions f(TZ). For 
example, let us consider the same matrix elements ()39j) but in THO representation: 



U q 



dz / rdr U q (z, r) 



f\K) df{K) 
-R? dK 



< 



zf(K) 
K 



< 



2 fr 2 nn) 
n 2 



(42) 



Introducing new dimensionless variables 



£ 



b z n 



r 2 pin) 

b 2 n 2 



(43) 



for which we have 



d£ di] 



bM 



f\K) df(K) 

n 2 dn 



rdr dz. 



(44) 



the matrix elements have the form of integrals, which are exactly identical to those in the HO 
basis (J4*nj) . after changing the function U q (£,T]) to 



U 9 fov) = kU q (tb,-fa,Jrjb ± jfc) 



(45) 



The calculation of matrix elements corresponding to derivative terms in the Hamiltonian (|17|) 
can be performed in an analogous way, after the derivatives of the Jacobian, 7^2 J% ; axe 
taken into account. 
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3.7 Calculation of Local Densities 

After diagonalizing the HFB equation (J3*7j) . local densities are calculated as 

p(nx,rV) = Y. P*P $*(ra)$^(rV) , 
p(ra,rV) = E^ ^(rtr)^(rV) , (46) 

a/3 

where <J> Q (rcr) denotes the HO or THO basis wave functions, and the matrix elements of mean- 
field and pairing density matrices read 

Pa/3 = Y^ V »kVf3k, Pap = ~ ^2 V ak U Pk ■ (47) 

k k 

The HFB calculations for zero-range pairing interaction give divergent energies when in- 
creasing the number of quasiparticle states in the sums of Eq. (J4*7)l (see discussion in Ref. [3]). 
Therefore, they invariably require a truncation of quasiparticle basis by defining a cut-off quasi- 
particle energy and including all quasiparticle states only up to this value. 

The choice of an appropriate cut-off procedure has been discussed in j2]. After each iteration, 
performed with a given Fermi energy A, one calculates an equivalent spectrum e k and pairing 
gaps A k : 

e k = (1 - 2N k )E k , 

, (48) 

A k = 2E k ^N k (l-N k ), 

where N k denotes the norm (|23|) of the lower HFB wave function. Using this spectrum and 
pairing gaps, the Fermi energy is readjusted to obtain the correct value of particle number, and 
this new value is used in the next HFB iteration. 

Due to the similarity between the equivalent spectrum e k and the single-particle energies, 
one can take into account only those quasiparticle states for which 

where e max >0 is a parameter defining the amount of the positive-energy phase space taken into 
account. Since all hole-like quasiparticle states, N k <l/2, have negative values of e k , condition 
(149)1 guarantees that they are all taken into account. In this way, a global cut-off prescription 
is defined which fulfills the requirement of taking into account the positive-energy phase space 
as well as all quasiparticle states up to the highest hole-like quasiparticle energy. In the code, 
a default value of e max =60 MeV is used. 

3.8 Coulomb Interaction 

In the case of proton states, one has to add to the central potential the direct Coulomb field 

KfW=e 2 |dV^l, (50) 

as well as the exchange Coulomb field, which in the present implementation is treated within 
the Slater approximation: 

^(r) = -e 2 (i) 1/3 pi /3 (r). (51) 
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The integrand in the direct term (|5U|) has a logarithmic singularity at the point r=r'. A way 
to bypass this difficulty is to use the Vautherin prescription [3Uj, i.e., to employ the identity 

A r /|r-r'| =2/|r-r'|, (52) 

and then integrate by parts the integral in Eq. (|5U|) . As a results, one obtains a singularity- free 
expression 

V d c (r) = ^Jdh>\r-r'\A T ,p p (r'). (53) 

In cylindrical coordinates, after integrating over the azimuthal angle ip, one finds 

/<» poo / ^ rr t \ 

rdr J dz y/d(r,z) E I —— \ Ap p (r, z), (54) 

where d(r, z) = \{z — z') 2 + (r + r') 2 ] and E(x) is the complete elliptic integral of the second 
kind that can be approximated by a standard polynomial formula [3Tj . 

Equivalently, one can use the prescription developed originally for calculations with the 
finite-range (Gogny) force [3]. It consists of expressing the Coulomb force as a sum of Gaussians: 

(55) 
(56) 





1 


9 f°° All. (-r') 2 








|r — r'| 


V^ Jo V 2 l 


which gives 


V d C (r) = 


^f>>> 


where the integral 


W = 


f _ (r-r') 2 




/ dVe ^ p(r') 



(57) 

can be easily calculated in cylindrical coordinates. After integrating over the azimuthal angle 
if, one finds 

I^r', z') = 2tt / rdr dz e ? J f — r J Pp (r, z), (58) 

where Iq(x) is the Bessel function that can also be approximated by a standard polynomial 
formula [3*Tj . 

In order to perform the remaining one-dimensional integration in Eq. ()56j) . the variable /i 
is changed to 

£ = 6/V^+^, (59) 

where 6 is the largest of the two HO lengths 6 Z and 6^. This change of variable is very convenient, 
since then the range of integration becomes [0, 1]. The integral (|5l)j) is accurately computed by 
using a 30-point Gauss- Legendre quadrature with respect to £. 

We have tested the precision of both prescriptions, Eqs. (|5*3~j) and (J56|) . and checked that 
the second one gives better results within the adopted numbers of Gauss-Hermite and Gauss- 
Laguerre points that are used for calculating proton densities. Therefore, in the code HFBTHO 
(vl.66p) this second prescription is used, while the first one remains in the code, but is inactive. 
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3.9 Lipkin-Nogami Method 

The LN method constitutes an efficient method for approximately restoring the particle num- 
bers before variation [33]. With only a slight modification of the HFB procedure outlined above, 
it is possible to obtain a very good approximation for the optimal HFB state, on which exact 
particle number projection then has to be performed [34"] l35j. 

In more detail, the LN method is implemented by performing the HFB calculations with an 
additional term included in the HF Hamiltonian, 

ti = h-2\ 2 (l-2p), (60) 

and by iteratively calculating the parameter A2 (separately for neutrons and protons) so as to 
properly describe the curvature of the total energy as a function of particle number. For an 
arbitrary two-body interaction V, X 2 can be calculated from the particle-number dispersion 
according to [33], 

(0|?|4 >( 4|*»|0) 

(0|iV 2 |4)(4|iV 2 |0> 

where |0) is the quasiparticle vacuum, N is the particle number operator, and |4)(4| is the 
projection operator onto the 4-quasiparticle space. On evaluating all required matrix elements, 
one obtains J3B] 

4Trr'p(l-p)+4TrA'(l-p)K 

2 8 [Trp(l - p)f - 16Trp 2 (l - p) 2 ' l ' 

where the potentials 

r L' = J2pp> V af3a'l3'{p(l - P))/3'I3, 

(63) 

can be calculated in a full analogy to T and A by replacing p and k by p(l — p) and pn, 

respectively. In the case of the seniority-pairing interaction with strength G, Eq. (J62J1 simplifies 

to 

_G Tt(1-p)kTt P k-2Tt(1-p) 2 P 2 

2 4 [Trp(l - p)f - 2 Trp 2 (l - p) 2 ' l ' 

An explicit calculation of A2 from Eq. (J62|) requires calculating new sets of fields (|63|) . which 
is rather cumbersome. However, we have found |2S] that Eq. (|H2*|) can be well approximated by 
the seniority-pairing expression (|6~4*j) with the effective strength 

A 2 
G = G cS = --r- (65) 

-C'pair 



determined from the pairing energy 
and the average pairing gap 



£pair = -|TrAK (66) 



^^- <-) 



Such a procedure is implemented in the code HFBTHO (vl.66p). 
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3.10 Particle-Number Projection After Variation 

Introducing the particle-number projection operator for N particles, 

pN = _L f d<f> e ^- N \ (68) 



where N is the number operator, the average HFB energy of the particle-number projected 
state can be expressed as an integral over the gauge angle <fi of the Hamiltonian matrix elements 
between states with different gauge angles (HZl EH] • In particular, for the Skyrme-HFB method 
implemented here, the particle-number projected energy can be written as [2H 123] 

(§\HP N \§) f f 

E N [P, P] = ($ | P V| $) = J <W Vtt) J ^ W(r, > ( 69 ) 

where the gauge-angle dependent energy density Tt(r, <p) is derived from the unprojected energy 
density TC(r) (JTU]) by simply substituting the particle and pairing local densities p(r), p(r), r(r), 
and Jjj(r) by their gauge-angle dependent counterparts p(r, (f>), p(r, (f>), r(r, 4>), and J»j(r, </>), re- 
spectively. The latter densities are calculated from the gauge-angle dependent density matrices 
as 

p(nr,rV,0) = ^p act >(0) $*,(r', a')$ a (r, a), 

aa' 

(70) 
p(r<r,r'a',4>) = ^~p aa >{4>) $* a ,{r',a')® a (r,a), 

aa' 

where the gauge-angle dependent matrix elements read 

Pa'a{4>) = y^Ca/3(0)P/3«' i 

P 

(71) 

Pa'a((p) = e" 40 22 C af3((fi)P(3a>, 
9 

and depend on the unprojected matrix elements (pT7j) and on the gauge-angle dependent matrix 

C(<f)) = e 2i<t> [1 + p(e 2i(t> - 1)] _1 . (72) 

Function y{<p) appearing in Eq. (JoU|) is defined as 

x(d>) r ,,. , e-^ N det(e^I) , N 

y U) = - vy ; - for x(0) = f- v ; , 73 

yKV> ftyxtf) W 27r ^detC(0) l J 

where J is the unit matrix. 

Since the gauge-angle dependent matrices (|70]) and (J7TJ) are all diagonal in the same canoni- 
cal basis that diagonalizes the unprojected density matrices (|37]> . all calculations are very much 
simplified when they are performed in the canonical basis. In particular, in the canonical basis 
the matrices (|7T]) read 

e 2t ^v 2 e^u v 

PM = 2 i 2zA 2 aIld PM = 2 | ft,f 2 . ( 74 ) 

iff, + e 2l fvf, u,, + e ltc fvi 

p A 1 f* P 
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while the function x{4>) can be calculated as 

*(0) = e-^nK + e2 ^)' ( 75 ) 

where v^ and u^ (v 2 + u 2 l = 1) are the usual canonical basis occupation amplitudes. 

All the above expressions apply to independently restoring the proton and neutron numbers, 
so, in practice, integrations over two gauge angles have to be simultaneously implemented. In 
practice, these integrations are carried out by using a simple discretization method, which 
amounts to approximating the projection operator (|55|) by a double sum J3H], i-e., 



where 



l-\ 1 L-\ 

pNz = _L V^ e i<p„(N~N)±_ y^ e i<t> P (z-z)^ ^ 

i n =a i P =o 



TV 

-lq, q = n,p. (77) 



Yq L 
Usually no more than L = 9 points are required for a precise particle number restoration. 

3.11 Constraints 

In the code HFBTHO (vl.66p), the HFB energy (JHJ) can be minimized under the constraint of a 
fixed quadrupole moment. This option should be used if one is interested in the potential energy 
surface of a nucleus along the quadrupole collective coordinate. The quadrupole constraint is 
assumed in the standard quadratic form [10] : 

£ q = C q ((Q)-q) 2 , (78) 

where (Q) is the average value of the mass-quadrupole-moment operator, 

Q = 2z 2 - r 2 , (79) 

Q is the constraint value of the quadrupole moment, and Cq is the stiffness constant. 

4 Program HFBTHO (vl.66p) 

The code HFBTHO (vl.66p) is written in Fortran 95 with MODULE definitions that specify all 
common arrays and variables for other subroutines by using the USE statements. Integer and 
real types of variables are automatically detected for the particular computer through the KIND 
statements. The code is entirely portable. It contains all initial data and no references to 
external subroutines or libraries are made. 

The code requires one input data file (tho.dat). Optionally, in case one wants to restart 
calculations from a previous run, two more files, dnnn_zzz.hel and/or dnnn_zzz.tel, are 
required as described below. Also optionally, if one wants to run the code for user-defined 
Skyrme-force parameters, file forces.dat is required. 

The results are printed on the standard output and also recorded in the file thoout.dat. 
The main results are also recorded in the files hodef.dat (HO basis) and thodef.dat (THO 
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basis), where one line is written for every nucleus calculated, producing a concise table of results 
suitable for further analyses. Files hodef . dat and thodef . dat are also used when restarting 
the given calculation after an abnormal termination (CPU time limit or system crash). Namely, 
before performing a given run, the code always checks if the line corresponding to this run is 
present or not in the file hodef.dat or thodef.dat. If this is the case, the code does not 
repeat the calculation for the given run, and only the runs which have not been completed are 
executed. Due to this implementation, if the user wishes to rerun the same input data file, files 
hodef . dat and thodef . dat have to be first removed from the current directory. 

4.1 General Structure of the Code 

The code runs, in sequence, the set of main subroutines listed in Table Q] If multiple runs are 
requested in a single input data file, the code always repeats the whole sequence of calls from 
the beginning to end, including an initialization of all variables and data. 

4.2 Input Data File 

Input data are read from file tho.dat, which is shown in Tabled The file consists of the first 
line, which contains only two numbers, below referred to as II and 12, followed by a sequence of 
identical lines, each of them defining one specific run of the code. All numbers containing a dot 
are type real, and those without a dot are type integer. In quotations there are four-character 
strings giving acronyms of the Skyrme forces. The code uses free format, so at least one space 
is needed in order to separate the input numbers. 
The code has three main regimes: 

(i) nucleus-after-nucleus, defined by I1<0, 
(ii) file-after-file, defined by 11=0, 

(iii) chain-after-chain, defined by II >0. 

In the nucleus-after-nucleus regime, the code ignores the values of |I1| and 12, and then 
performs one run for each line of the input data file that follows the first line. This is the 
simplest and most often used regime, illustrated by the example given in Table [T] 

In the file-after-file regime, the code ignores the value of 12, and then reads the second line 
of the input data file, from where it takes all fields except from the values of ININ, N, and Z. 
Then it performs one run for each dnnn_zzz file found in the current directory. Files dnnn_zzz 
contain results of previous runs and are described below. 

In the chain-after-chain regime, the code reads the second line of the input data file, from 
where it takes all fields except from the values of N, and Z. Then it performs one run for each 
nucleus in the chain of isotones or isotopes located between the bottom of the stability valley 
and the drip line. The bottom of the stability valley is parametrically defined as 

f(N, Z) = N - Z - 0.006(iV + Z) 5/3 = 0. (80) 

• If I2>0, the code calculates the chain of isotopes for the proton number Z=ll, starting 
with the lowest even neutron number N satisfying f(N,Z)>0, and then step-by-step 
increasing the number of neutrons by two. Calculations continue until the neutron drip 
line is reached, and then the program stops. 
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Table 1: List of main subroutines constituting the code HFBTHO (vl.66p) 



Subroutine 



Task 



DEFAULT 
READ INPUT 
PREPARER 

BASEO 

THOALLOC 
BASE 

GAUPOL 

INOUT 

ITER 



F01234 
ITER 

RESU 
INOUT 



Initializes all variables (initially, or after the previous run). 

Reads parameters from the input data file tho.dat. 

Initializes variables according to the user's request defined in the input 

data file. 

Determines the HO configurational space and dimensions of allocat- 

able arrays. 

Allocates memory required for the given run of the code. 

Calculates and stores properties of the configurational space and all 

associated quantum numbers. 

Calculates and stores the HO basis wave functions. 

Sets or reads (optional) initial densities, fields, and matrix elements. 

Main iteration loop for the HFB+HO calculation, which is repeated 

until convergence is met. It includes the following subroutines: 

DENSIT Calculates densities in coordinate space. 

FIELD Calculates mean fields in coordinate space. 

GAMDEL Calculates the particle-hole and pairing Hamiltonian 

matrices. 
EXPECT Calculates average values of observables. 
HFBDIAG Diagonalizes the HFB equation. 

After the HFB+HO solution is found, calculates the THO basis wave 

functions, which replace the HO ones. 

Main iteration loop for the HFB+THO calculation, which is repeated 

until convergence is met. The same subroutine and sequence of calls 

is used as above. 

Calculates all required physical characteristics and canonical basis 

properties, and performs the particle number projection. 

Records the final densities, fields, and matrix elements for feature use 

(optional). 



Table 2: Input data file tho.dat. 



(a) 


(b) 


(c) 


(d) (e) 


(f) 


(g) 


(h) 


(i) 


0) 


(k) 


(1) (m) (n) 


(o) 


(P) 


(q) 


W (s) 


-20 


20 




























20 


-2. 


0. 


-1 300 


1 


70 


50 


'SLY4' 


-1 


1 


0.26 0.5 


9 





2 


2 0.0001 


20 


-2. 


0. 


-1 300 


1 


72 


50 


'SLY4' 


-1 


1 


0.26 0.5 


9 





2 


-2 0.0001 


20 


-2. 


0. 


-1 300 


1 


74 


50 


'SLY4' 


-1 


1 


0.26 0.5 


9 





-2 


2 0.0001 


20 


-2. 


0. 


-1 300 


1 


76 


50 


'SLY4' 


-1 


1 


0.26 0.5 


9 





-2 


-2 0.0001 


-14 


-2. 


0. 


-1 300 


1 


78 


50 


'SKP ' 


-1 


1 


0.26 0.5 


9 





4 


4 0.0001 





-2. 


0. 


-1 300 


1 


70 


50 


'SLY4' 


-1 


1 


0.26 0.5 


9 








0.0001 
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• If I2<0, the code calculates the chain of isotones for the neutron number N=ll, starting 
with the lowest even proton number Z satisfying f(N,Z)<0, and then step- by-step in- 
creasing the number of protons by two. Calculations continue until the proton drip line 
is reached, and then the program stops. 

All lines of the input data file, after the first line, contain 19 fields each. Below we denote 
these fields by letters (a) - (s), as shown in the header of Tabled The description of the fields 
is as follows: 

• (a) Number of oscillator shells N s h'. 

— If N s h > 0, the code prints intermediate results at every iteration. 

— If N s h < 0, the code prints results at the first and last iterations only, and the module 
of the input value is used for N s h- 

— If N s h = 0, the code stops. This value is used to indicate the end of the input data 
file. 

For N s h > 14, the code always begins with a short, 20-iteration run using N s h = 14, and 
the resulting fields then serve as a starting point for the calculation with the requested 
value of N sh . For the THO-basis calculations, use of N sh < 14 is not recommended, 
because precision of the HO density profile can be insufficient for a reliable determination 
of the LST function. 



(b) Oscillator basis parameter & =V"T+^l : 

— If &o > 0, the code uses this given value of bo- 

— If b < 0, the code uses the default value of b Q = ^2(h 2 /2m)/(41fA- 1 / 3 ) for f=1.2. 

(c) Deformation /3 of the HO basis. The value of /3q defines the HO oscillator lengths 
through 6j_ = b q~ 1 / 6 , b z = &o<? 1//3 > an d q = exp (3^5/(16ir)P ). In particular, the value 
of Po = corresponds to the spherical HO basis. 

(d) The THO basis control parameter ILST: 

— If ILST= 0, the code performs the HO basis calculation only. If ININ< 0, the file 
dnnn_zzz.hel is used as the starting point. If MAXI>0, at the end of the given run 
file dnnn_zzz.hel is stored. 

— If ILST= — 1, the code performs the HO basis calculation followed by the THO 
basis calculation. If ININ< 0, the file dnnn_zzz.hel is used as the starting point. 
If MAXI>0, at the end of the given run files dnnn_zzz.hel and dnnn_zzz.tel are 
stored. 

— If ILST= 1, the code performs the THO basis calculation only. File dnnn_zzz.tel 
must exist and is used as the starting point (only ININ<0 is allowed). If MAXI>0, at 
the end of the given run file dnnn_zzz .tel is stored. 
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(e) Maximal number of iterations MAXI. If the negative number is read, the absolute value 
is used. 

— If MAXI>0, at the end of the given run files dnnn_zzz.hel and/or dnnn_zzz.tel are 
stored, 

— If MAXKO, files dnnn_zzz.hel and dnnn_zzz.tel are not stored, and the module of 
the input value is used for MAXI. 

(f) The starting-point control parameter ININ: 

— If ININ= 1, the code starts from a default spherical field predefined within the code, 

— If ININ= 2, the code starts from a default prolate field predefined within the code, 

— If ININ= 3, the code starts from a default oblate field predefined within the code, 

— If ININ= —1, the code starts from file snnn_zzz.hel or snnn_zzz.tel, 

— If ININ= —2, the code starts from file pnnn_zzz.hel or pnnn_zzz.tel, 

— If ININ= —3, the code starts from file onnn_zzz.hel or onnn_zzz.tel. 

(g) Number of neutrons N. 

(h) Number of protons Z. 

(i) Skyrme force character*4 acronym, e.g., 'SIII', 'SKP ', 'SLY4', or 'SKM*'. If value 
'READ' is read, the code reads the Skyrme force parameters from file forces.dat. An 
example of the file forces.dat is presented in Table 01 

(j) The Lipkin-Nogami control parameter KINDHFB: 

— If KINDHFB= 1, Lipkin-Nogami correction not included, 

— If KINDHFB=— 1, Lipkin-Nogami correction included. 

(k) The pairing- force control parameter IPPFORCE: 

— If IPPF0RCE= 0, No pairing correlations (Hartree-Fock calculation), 

— If IPPF0RCE= 1, Calculation for the density-dependent delta pairing force, 

— If IPPF0RCE= 2, Calculation for the density-independent delta pairing force. 

(1) The quadrupole-constraint control parameter ICSTR. If ICSTR=0, the quadrupole con- 
straint is not included, and the next two fields (m) and (n) are not used. If ICSTR=1, 
then: 

— (m) Constrained value of the quadrupole deformation (3. The value of j3 defines the 
constrained quadrupole moment Q in Eq. (|78JI through: Q = J ~ < r 2 > J3. To be 
done: what is < r 2 >. 

— (n) Parameter r\ defining the stiffness Cq of the quadratic quadrupole constraint 
constant by C Q = J7(4LA -1 / 3 )/(8A&g < r 2 >). 
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Table 3: User-defined parameters of the Skyrme force, as given in the file forces.dat. 



Value 


Description 


'SLY4' 



-0.2488913d+04 


Skyrme-force acronym 

Tensor term (O-excluded, 1-included) 

to 


0.4868180d+03 


h 


-0.5463950d+03 


t 2 


0.1377700d+05 


h 


0.8340000d0 


X 


-0.3440000d0 


Xi 


-l.OOOOOOOdO 


X2 


1.3540000d0 


%3 


0.1230000d+03 


w 


6.0d0 


1/a 


20.735530d0 


h 2 /2m 


0.160d0 

l.OdO 

60.0d0 

0.5d0 


Po (saturation density for pairing) 
7 (power of density for pairing) 
e ma x (pairing cut-off energy) 
V\ (0-volume, 1-surface, 0.5-mixed) 


-244.7200d0 


Vq (pairing strength) 



• (o) The number of gauge-angle points L used for the particle number projection. Note 
that the code always performs the PNP, even if pairing correlations are not included. 

• (p) The particle-number-shift control parameter ISHIFT. If, ISHIFT=0, the particle- 
number-projection is performed on N and Z, and the next two fields (q) and (r) are 
not used. If ISHIFT=1, then: 

— (q) Neutron- number shift KDN, i.e., the projection is performed on neutron number 

7V+KDN, 

— (r) Proton- number shift KDZ, i.e., the projection is performed on proton number 
Z+KDZ. 

• (s) Requested precision of convergence SI (in MeV). Iterations stop when changes of 
all mean-field and pairing matrix elements between two consecutive iterations become 
smaller than the value of SI. Recommended value is 0.0001. 

After the solution is found, and if MAXI>0, the code stores files dnnn_zzz.hel (if the HO- 
basis run has been performed) and/or dnnn_zzz . tel (if the THO-basis run has been performed). 
Names of these files are automatically constructed based on the input-data parameters ININ, 
N, and Z, namely: 



V, 'p', or V, 



for |ININ| = 1, 2, or 3, respectively, 



21 



• nnn = three-digit value of N with leading zeros included, 

• zzz = three-digit value of Z with leading zeros included. 

These files can be used in a later run to restart calculations from previously found solutions. 
For example, file s070_050.tel contains results of the THO-basis calculation for 120 Sn, which 
has been obtained by starting from a spherical field. Note that the name of the file reflects the 
starting deformation only, while it may, in fact, contain results for another deformation that 
has been obtained during the iteration. 

4.3 Output Files 

The results are printed on the standard output file. Each run produces a separate part of the 
output file; also the HO run preceding a THO run produces one such part. Below we briefly 
describe different sections of the output file. 

• Header. Contains the version number of the code, date and time of execution, name of 
the element, and its particle, neutron, and proton numbers. 

• Input data. Contains a short summary of the input data for the requested run. 

• Force. Lists the acronym and parameters of the Skyrme force, as well as parameters of 
the pairing force. 

• Numerical. Contains some information on numerical parameters and options used for the 
given run. 

• Regime. Gives the regime in which the code is run. 

• Iterations. Shows brief information about iterations performed. One line of the output 
file per each iteration is printed and contains the following columns: 

— Iteration number i. 

— Accuracy si. 

— Current mixing parameter between the previous and current fields mix. 

— Quadrupole deformation beta, /3 = a/|'<^>' ^ or Q gi ven m Eq. (|75j). 

— Total energy Etot. 

— Particle number A. 

— Neutron rms radius rn. 

— Proton rms radius rp. 

— Neutron pairing energy En. 

— Neutron pairing gap Dn. 

— Proton pairing energy Ep. 

— Proton pairing gap Dp. 
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— Neutron Fermi energy Ln. 

— Proton Fermi energy Lp. 

• Files. Contains information on the dnnn_zzz.hel or dnnn_zzz.tel file written. 

• Observables. Lists values of various observables calculated for the HFB state without 
PNP and with the Lipkin-Nogami corrections, and then those calculated for the PNP 
HFB state. 

The same information, plus more results on the quasiparticle and canonical states, is also 
stored in the file thoout .dat. However, this file is rewound after each run, so it contains results 
of only the last run executed for the given input data file. 

Files hodef . dat and thodef . dat contain synthetic results of all runs, printed in the form of 
a single line per each performed run. If the given run performs only an HO-basis calculation, or 
only a THO-basis calculation, then only an entry in file hodef . dat or thodef . dat is produced, 
respectively. On the other hand, runs that perform both HO and THO calculations produce 
entries in both these files. Lines in the files hodef.dat and thodef.dat contain 105 columns 
each, and each column is described by a name printed in the first header line. The names 
are self-explanatory, and most often they correspond to the names used in the present write- 
up. Names preceded by U: pertain to results obtained for the HFB states before PNP, while 
those beginning with L pertain to the results containing the Lipkin-Nogami corrections. Names 
ending with t, n, or p give total, neutron, or proton observables, respectively. 

5 Conclusions 

The code HFBTHO (vl.66p) is a tool of choice for self-consistent calculations for a large number 
of even-even nuclei. Several examples of deformed HFBTHO calculations, recently implemented 
on parallel computers, are given in Ref. """"""J- By creating a simple load-balancing routine 
that allows one to scale the problem to 200 processors, it was possible to calculate the entire 
deformed even-even nuclear mass table in a single 24 wall-clock hour run (or approximately 
4,800 processor hours). 

The crucial input for such calculations, which determines the quality of results, is the nu- 
clear energy density functional. The development of the "universal" nuclear energy density 
functional still remains one of the major challenges for nuclear theory. While self-consistent 
HFB methods have already achieved a level of sophistication and precision which allows analy- 
ses of experimental data for a wide range of properties and for arbitrarily heavy nuclei (see, e.g., 
Refs. |4*T| 14*2*1 14*5] for deformed HFB mass table), much work remains to be done. Developing 
a universal nuclear density functional will require a better understanding of the density de- 
pendence, isospin effects, and pairing, as well as an improved treatment of symmetry-breaking 
effects and many-body correlations. 

In addition to systematic improvements of the nuclear energy density functional, there are 
several anticipated extensions of HFBTHO itself. The future enhancements to HFBTHO will 
include the implementation of the full particle-number projection before variation, extension of 
code to odd particle numbers, implementation of non-standard spin-orbit term and two-body 
center-of-mass correction, and evaluation of dynamical corrections representing correlations 
beyond the mean field. 
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